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We apply the event-chain Monte Carlo algorithm to classical continuum spin models on a lattice 
and clarify the condition for its validity. In the two-dimensional XY model, it outperforms the 
local Monte Carlo algorithm by two orders of magnitude, although it remains slower than the Wolff 
cluster algorithm. In the three-dimensional XY spin glass model at low temperature, the event-chain 
algorithm is far superior to the other algorithms. 
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INTRODUCTION 


Classical and quantum spin models are of fundamen¬ 
tal interest in statistical and condensed-matter physics. 
Spin models are also a crucial test bed for computational 
algorithms. 

An important representative is the model of con¬ 
tinuous two-dimensional classical spins of fixed length 
(rotators) on a two-dimensional lattice. Thirty years 
ago, the existence and nature of the phase transi¬ 
tion in this two-dimensional XY model were highly 
controversial pLj. The substitution of the traditional lo¬ 
cal Monte Carlo (LMC) algorithm^ by Wolff’s spin flip 
cluster (SFC) algorithm |3] then quickly allowed to clarify 
that this model indeed undergoes a Kosterlitz-Thouless 
transition’5], whose temperature is now known to five 
significant digits m d. SFC has played a decisive role in 
understanding the physics of the XY model [8UT0]. and 
in arriving at its detailed quantitative description. 

SFC and its variants can be implemented for a wide 
range of models, but they are efficient only in a few of 
them. Particularly frustrating is the case of the three- 
dimensional XY spin glass model, where the algorithm 
loses all its power jTTl [12] . For this much studied spin 
glass model, our understanding today resembles the one 
of the XY model before the revolution triggered by the 
cluster algorithms. Clearly, there still is a great need for 
more powerful algorithms for classical and quantum spin 
models. 

Today’s Markov-chain Monte Carlo algorithms gen¬ 
erally follow the conventional paradigm based on three 
principles: 1/ Each move represents a finite change of 
the configuration. It is independent of the previous move, 
and depends only on the configuration itself. 2/ The al¬ 
gorithm satisfies the detailed-balance condition. 3/ The 
decision whether a proposed move is accepted is based 
on the change in energy, using the Metropolis acceptance 
rule or the heat-bath conditional [131. 

In the present work, we show that the novel event- 
chain Monte Carlo (ECMC) paradigm [T4lfl6] . that has 
already been very successful in particle systems mm , 


can also be applied to the XY model and the XY spin 
glass model. The paradigm breaks all three principles 
of the conventional Markov-chain scheme: Moves are in¬ 
finitesimal rather than finite, although an event-driven 
scheme allows to recover finite displacements [16] . In one¬ 
dimensional systems, the moves do not change with time. 
In multidimensional systems, moves persist on long time 
scales. This is achieved within the Markov-chain scheme 
through additional “lifting” variables [15] [21]. In addi¬ 
tion, ECMC violates detailed balance and only satisfies 
the weaker global balance condition (cf. [22H26] ). Fi¬ 
nally, the decision on future moves is based on the change 
in pair energies, rather than the change in total energy. 
This is achieved by replacing the standard Metropolis al¬ 
gorithm by its recently introduced factorized variant fT5|. 

For the two-dimensional XY model at the critical 
point, we find that ECMC is about 100 times faster than 
LMC, although the presence of a slow time scale in auto¬ 
correlation functions makes that it is not as fast as SFC. 
In the low-temperature phase of the three-dimensional 
XY spin glass model, where SFC is known to be ineffi¬ 
cient, ECMC clearly outperforms LMC. 


FROM LOCAL MONTE CARLO TO THE 
“EVENT-CHAIN” ALGORITHM 
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FIG. 1. LMC move for the one-dimensional XY model. Up¬ 
per panel : Configuration at time t and proposed displace¬ 
ment A<f of a randomly chosen spin, corresponding to an en¬ 
ergy change A E. Lower panel : Possible configurations at 
time t+1: The proposed move is accepted with probability 
min(l, exp(— /3AE) (left) and rejected otherwise (right). 


In the two-dimensional ferromagnetic XY model of 
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spins Sk = (5^,5^) = (cos </>&, sin cj)^) on a lattice with 
sites i = 1,..., TV, and with an energy 

E = ^ ^ JijSi • iSj = ^ ^ [ Jij cos(02 0j)L (1) 

(hi) (hi) ^ ^ 

^13 

the coupling constants Jij are all equal to one. The sum 
(i,j) goes over nearest neighbors on the lattice. We refer 
to the Eij as “pair energies”. The XY model on a two- 
dimensional square lattice undergoes a phase transition 
at inverse temperature f3 = 1.1199, see ref. [6|. 

In LMC, one proposes at each time step t a finite move 
from a configuration a to a configuration b (a rotation by 
a finite angle Acj) of a spin fc), as sketched in Fig. |T] 
To satisfy detailed balance [13], k is randomly chosen at 
each time step, and Ac/) is sampled from a symmetric 
distribution around zero, so that Acj) arises with the same 
probability as — Acj). The proposed move corresponds to 
an energy change A E = E^ — E a in Eq. [lj and it is 
accepted with probability 

^ = min(l,exp(-/3A£)). (2) 

The exponential in this equation corresponds to the ratio 
^b/^a of the Boltzmann weights of the configurations. 

Practically, the move is accepted, and the configuration 
updated to 6, if a uniform random number between 0 
and 1 satisfies ran(0,1) < (see [13j). Otherwise, the 
configuration at time t +1 is the same as the one at time 
£, namely a. 

The recently introduced factorized algorithm [151 also 
satisfies the detailed-balance condition. In this method, 
the energy-based Metropolis acceptance probability is re¬ 
placed by a factorized form which separately depends on 
the pair-energy changes: 

Pa“c = II Pacc = II min(l,exp(-/3A£fe/)). (3) 

(k,l) (k,l) 

The proposed move a b is accepted with this probabil¬ 
ity. The factorized algorithm always has a smaller accep¬ 
tance rate than the conventional one, p^c < P^ (this 
will however turn out not to be a problem in ECMC). To 
implement Eq. [3j one might use a single random num¬ 
ber and accept the move if ran(0,1) < p \We rather 
accept the move if several independent random numbers 
satisfy ran/ c /(0,1) < p^ l cc for all pairs k,l. In other words, 
a move is accepted only if it is pair-accepted by all pairs 
k,l. This consensus rule is illustrated in Fig. [2] We 
note that the factorization in Eq. [3] relies on the pos¬ 
sibility to cut the hamiltonian into independent pieces. 
The factorization may also be used to separate different 
components of the inter-particle potential, as for example 
the 1/r 6 and 1/r 12 pieces in the Lennard-Jones potential 

nsiiis]. 

The ECMC combines the factorized Metropolis prob¬ 
ability with the “lifting” concept of Diaconis et al. m 
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FIG. 2. Factorized Metropolis move. Upper panel : Configu¬ 
ration at time t and proposed displacement Acj) of a randomly 
chosen spin k. Middle panel : Factorization into pairs (j, k) 
and (k, l). In the factor ( j , k), the move is pair-accepted with 
probability min(1, exp (—/3AEjk), etc. Lower panel: Possible 
configurations at time t + 1: The proposed move is either ac¬ 
cepted by consensus (i.e. independently by all pairs) or else 
rejected. 


and with the idea of infinitesimal displacements [13]. The 
term “lifting” refers to the extension of the physical con¬ 
figuration by an additional variable that fixes the pro¬ 
posed move. Written as it singles out the spin k as 
the only one that can move, as fa fa Y Acj) (see Fig.[3|. 
If the move is accepted, the lifting variable for the next 
time step t + 1 is again & . If the physical move is re¬ 
jected, a lifting move takes place and the lifting variable 
is passed on to the spin l of the pair that rejected the 
move, and the physical configuration is unchanged. In 
both cases, the value of Acj) is used again. Note that 
for infinitesimal A0, the acceptance probabilities of the 
physical moves approach one and the rejection probabil¬ 
ities approach zero. Multiple rejections are totally sup¬ 
pressed, and the choice of ^ is unique [15]. At each time 
step, either a lifting move or a physical move takes place, 
and ECMC is thus formally reject ion-free. 

ECMC satisfies the global balance condition in the 
XY model, as we now show: For simplicity, we consider 
only two spins and concentrate on a configuration d (see 
Fig.§. This configuration can only be reached through 
a lifting move from a or through a physical move from 
b. The global-balance condition p~3] states that the flow 
into configuration d must be equal to the flow out of it: 

7 r a p(a -A d) + 7 Tbp(b d) = 

v -V-' v -V-' 

V(a—>d) V(b—±d) 

7c d p(d -> /) + Tc d p(d -A a) . (4) 

V(d^f) V(d-)-a) 

Here, V(a d) represents the probability flow from a 
to d, etc. For ECMC, the probabilities p in Eq. [4] coin¬ 
cide with the acceptance probabilities: All configurations 
carry a lifting variable that specifies the spin that may 
move and the move itself, Acj). 

The statistical weight 7r a is trivially equal to ir d be- 
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FIG. 3. Lifting approach of ECMC. Physical moves b —>• d, 
d f and a —> c are by the same infinitesimal angle A(p in 
clockwise direction, all others are lifting moves that preserve 
the physical configuration. Note that 7 — tt c because of 

Eq.0 


cause they differ only by a lifting move. Furthermore, 
7r c equals 7p>, as the two configurations differ only by a 
global rotation. Writing A E = E^ — E d: we thus find 

V(b —> <f) = Tr b p { ^(b ~^d) = iT d Pa?c(d b) 

= ir d min(l, exp(-/3A£')). (5) 

Note in this equation that tt bPacci^ d) = KdP*£££{d 

b ), because the factorized transition probabilities satisfy 
detailed balance. Likewise, the change in energy in going 
from a c is also A E and p(a d) = 1 — p(a -Y c ). 
Therefore, the flow V(a -Y d) satisfies 

V(a -Y d) = 7r a (l — min(l,exp(— f3AE)) 

= 7r d (l - min(l, exp(-/JA2S)). (6) 

It follows that the flow into d, namely the sum of P(a -Y 
d ) and of P{b -Y d), equals 7r^. As for the flow out 
of d, it trivially equals n d because of the conservation of 
probabilities. It follows that the global balance of Eq.[4]is 
satisfied. The factorization property and the infinitesimal 
limit guarantee that the argument carries over to general 
N (see [H]). 

ECMC violates the detailed balance condition V(b -Y 
d) = V(d -Y b): A move d -Y b would be anti-clockwise, 
yet all moves within ECMC are, by the initial choice of 
A0, clockwise. Also, V(a -Y d) = 0, as E d > Ef and all 
physical moves from d to f are accepted. Furthermore, 
for ECMC to be valid, the pair energy must be symmetric 
(so that 7T5 = 7r c in Fig. [3|. Modified XY models, as 
described in ref. m, can also be treated, but more general 
pair energies require special considerations [28] . 

ECMC with infinitesimal moves requires a scaling of 
physical time: In one unit of time, as A 4 goes to zero, 
an infinite number of physical moves take place, but the 
number of lifting moves remains finite. In an event- 
driven approach [H, 16], the algorithmic complexity can 
be made to scale with the number of liftings: The lifting 
variable being set to the angle 4k now rotates clock¬ 
wise until the “event”, i.e. a lifting move, is produced 


through a rejection by a neighbor 1. The lifting variable 
is updated to 4i rotates clockwise, etc. Effectively, 
one undergoes an infinite number of Monte Carlo steps, 
giving a continuous trajectory. 

The angle 4k corresponding to the next event is easily 
sampled: We continue to consider a single pair (fc, l) of 
spins, with the lifting variable The i-th infinitesimal 
update of 4k is noted as the move i — 1 -Y i and the 
weight of the configuration (4i = 4k + zd0, 0z), 7q. The 
probability p event (0 -Y n) to accept n subsequent physical 
moves and then to reject the n Y 1st physical move is 


Pevent(0 ^ — Pacc(0 ^ 1) ' P&ccijl 1 ^ Tl) 

[1 — P&cc{n n + 1)] • (7) 

The j th term in this expression is min(l, nj/nj-i). Sup¬ 
posing for a moment that i\j is monotonously decreasing 
with j, this gives 


Pevent(0^n) = ^fl-^) 

TTO V Kn-l J 


— 1 dir 
7TO 94k 



(8) 


This probability is normalized, writing </> eve nt the value 
of 4k at which the event happens: 


1 

7T0 



dn 

d<pk 


dfieve nt 


—^event 


1 

7T0 



dir 


event 


= 1. 


( 9 ) 


This integral is sampled by m 


T^event — ran(0,7To) 
^’event/^’o = ran(0,1), 


( 10 ) 


which is equivalent to the following sampling of the en¬ 
ergy increase: 


A E(4 event ) = - [log ran(0,1)] f/3. (11) 


Sampling 7 r uniformly between 0 and the present value, 
7To (equivalently, A E from its exponential distribution) 
thus yields the event time, 4 event (see Fig. [4]). 

For a non-monotonous probability distribution, all 
negative energy increments correspond to an acceptance 
probability 1, and disappear from Eq. [ 7 ] The sampling 
of the energy increase in Eq. [TT] turns into the sampling 
of only the positive energy changes. As shown in Fig. [4] 
this can be expressed as a function E *, constructed only 
from the positive increments of the energy E [16] . 

For a system of more than one pair of spins, the event 
times 0event f° r each neighbor of the lifted spin k can be 
computed independently in view of the factorized prob¬ 
ability of Eq. [3] and k turns clockwise up to the earliest 
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^event = rail(0, 7T 0 ) 


<vent = ran(0,7rj) 




■^event [F§ bevent )] /ft 





FIG. 4. Event-driven implementation of ECMC for a pair 
of spins (k,l). From a starting point fa — fa of weight no 
and energy Eq, fa is updated by infinitesimal moves until 
4>k = favent- Left: Monotonously decreasing distribution 7 r: 
The lifting event is sampled as 7r eV ent = ran(0,7To). Right: 
General distribution 7r: E * vent — E*( 0) = [—log ran(0,1 )\//3. 


event (that involves, say, another spin l). The lifting 
variable is then set to 

It follows from Eq. [7] that all configurations encoun¬ 
tered between two events sample the Boltzmann distri¬ 
bution. Any uniform subset of these configurations can 
be used for averaging observables. A practical choice con¬ 
sists in outputting spin configurations at regular intervals 
independent of the occurrence of events. 

For the models considered here, we found that the 
efficiency was not increased by halting and restarting 
the simulation after fixed displacements. In contrast, 
switches between moves along the different coordinate 
axes assure ergodicity in multi-dimensional hamiltonians 
as they appear in particles systems [M , but also the re¬ 
lated Heisenberg model [29 . 


SIMULATIONS FOR THE TWO-DIMENSIONAL 
XY MODEL AT THE CRITICAL POINT 

In the two-dimensional XY model, we consider the 
susceptibility x 


v IIESfcll 2 
x N 


( 12 ) 


and estimate the convergence properties by the suscepti¬ 
bility autocorrelation function 


C x (t) = 


(xft / + ^)xft'))-(x ) 2 

<X 2 > - <x> 2 


(13) 


at the critical point (3 = 1.1199 (see 0). We suppose 
that x is a slow variable of this model. We measure time 


in sweeps: For ECMC, one sweep corresponds to ~ N 
lifting events while for LMC, one sweep corresponds to 
N attempted moves. For SFC, a sweep denotes ~ N 
spins added to clusters. The complexity of one sweep is 
O(N) in the three algorithms and the CPU times used 
per sweep are roughly comparable. 

In Fig. [5j we show the autocorrelation function for the 
XY model at its critical point, obtained from very long 
single runs of the algorithms. For LMC and SFC, the 
decay of the susceptibility autocorrelation function can 
be described by a single time scale, while for ECMC, it 
is well described by two time scales: 


"exp(—t/r LMC ) (LMC) 
exp(—t/r SFC ) (SFC) 

A 0 exp(—£/r FCMC ) + 
k Ai exp(-t/r 1 ECMC ) (ECMC) 


(14) 


For ECMC, this correlation function rapidly decays to 
^ 0.1 on a timescale r ECMC of about 5 sweeps. A 
slow mode r ECMC then sets in. It presents a 2 = 2 scal¬ 
ing ( r p CMC 00 L 2 , with N = L 2 ). As shown on the 
right panel of Fig. [5J r ECMC is an order of magnitude 
smaller than r LMC . Together with the initial rapid de¬ 
crease, this makes ECMC about one hundred times faster 
than LMC. However, its dynamical scaling exponent ap¬ 
pears to be z ~ 2, as for LMC. We notice that in particle 
systems, ECMC also shows initial ballistic behavior, but 
then crosses over into slower decay [501. 


THREE-DIMENSIONAL XY SPIN GLASS 
MODEL 


We now study ECMC for the three-dimensional XY 
spin glass model, where the nearest-neighbor coupling 
constants faj are drawn from a Gaussian normal distri¬ 
bution of zero mean and unit variance. The algorithm 
can be formulated as for the ferromagnetic model, and 
the spins continue to always turn clockwise. We will find 
evidence that the relaxation dynamics of ECMC differs 
from the one of LMC. Following m, we consider the 
chiral overlap between two independent systems, (1) and 
(2), with identical coupling constants 

a*) 

p =1 


with being the chirality at a plaquette p, perpen¬ 

dicular to the axis /i, defined as: 


(») _ 


1 

2\/2 


^2 s gn(J ii )sin(^ - 
(hj)£p 


(16) 


The sum ^ ep is taken over the four bonds encircling 
the plaquette p clockwise. By construction, p K is a sym¬ 
metric function about zero. As shown in Fig. [6j ECMC 
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FIG. 5. Autocorrelation function C x (t) for the two-dimensional XY model at the critical point [3 — 1. 119 9 for LMC (red, 


triangle ), ECMC (blue, circle ), and SFC (yellow, square). Exponential fits (black, dotted) are as in Eq. 14 Left: N = 32 


Middle: N = 128 2 . Right: Scaling of the autocorrelation time r with the system size. Both LMC (red, triangle) and the slow 
scale of ECMC (dark blue, circle) are compatible with a dynamical scaling exponent z ~ 2. Both the fast scale of ECMC (light 
blue, diamond) and SFC (yellow, square) are compatible with z ~ 0. Right Inset: Speedup of ECMC with respect to LMC vs. 
L. 


and LMC agree very well at high temperature. The au¬ 
tocorrelation function of the chiral overlap for LMC and 
ECMC, shown in Fig. |6j gives at high temperature a size- 
independent speedup by a factor ^ 5 of ECMC. 

The phase diagram of the three-dimensional XY spin 
glass model at low temperature (with the possible ex¬ 
istence of separate spin-glass and chiral-glass phases) is 
still being debated. We consider /3 = 3.636, which may 
be the locus of the spin glass transition [12], or below it, 
near the transition [3T1I32] . At this temperature, ECMC 
exhibits a striking advantage over LMC in one third of 
samples of size N = 6 3 , where it explores the configura¬ 
tion space more easily, without using parallel tempering 
[33] . A typical example of a symmetric chiral overlap dis¬ 
tribution profile after 10 6 sweeps (symmetric for ECMC, 
but not for LMC is shown in Fig. [ 7 ] together with the 
corresponding autocorrelation function. For larger sys¬ 
tems, the speedup of ECMC with respect to LMC seems 
to increase, but already for 10 3 systems, ECMC no longer 
equilibrates at (3 = 3.636. 


CONCLUSION 

In conclusion, we have applied in this work the re¬ 
cent event-chain algorithm to classical spin models, and 
obtained a considerable algorithmic speed-up with re¬ 
spect to the local Monte Carlo algorithm for the two- 
dimensional XY model at its critical point. The new 
method appears very general, as we also obtained clear 
acceleration for the three-dimensional XY spin glass 
model at low temperature. It will be interesting to see 
how well the event-chain algorithm couples with the tra¬ 
ditional acceleration methods, as for example the parallel 
tempering method, or the overrelaxation approaches that 



FIG. 6. Cumulative distribution of the chiral overlap p K 
for the three-dimensional XY spin glass model for N — 
4 3 ,6 3 ,8 3 ,10 3 at /3 = 1.5, in the high-temperature phase 
(single samples). Inset: Autocorrelation function C PK (t) for 
N = 6 3 from LMC (red, triangle) and ECMC (blue, circle). 


have been much used for spin glasses. 
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FIG. 7. Chiral overlap autocorrelation function from ECMC 
and LMC at /3 = 3.636 for a given sample at N = 6 3 . Inset: 
Distributions of p K after 10 6 sweeps for the two algorithms in 
the same sample. Note the nearly symmetric distribution for 
ECMC. 
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